
COMPILE_OPT IDL2, LOGICAL_PREDICATE

;; Load the PVEC structure!
conf = omni_load_conf()
IF ~exist(pvec) THEN restore,'./local/'+conf.survey+'_pvec.sav'

ndpdf = n_tags(pvec)-7
print,ndpdf

dnames = (tag_names(pvec))[1:ndpdf]

;;====================================
;; Table Header -- Simple Stats Table
openw,lun,'./grs_paper/dpdf_kdars.tex',/get_lun

printf,lun,'%%=================================%%'
printf,lun,'%%                                 %%'
printf,lun,'%%   Table: DPDF KDAR Table        %%'
printf,lun,'%%   Label: table:kdars            %%'
printf,lun,'%%                                 %%'
printf,lun,'%%=================================%%'
printf,lun,'\begin{deluxetable}{lcccc}'
printf,lun,'  \tablecolumns{5}'
printf,lun,'  \tablewidth{0pt}'
printf,lun,'  \tabletypesize{\footnotesize}'
printf,lun,'  \tablecaption{KDA Resolutions for BGPS V2 Sources\label{table:kdars}}'
printf,lun,'  \tablehead{'
printf,lun,'    \colhead{KDA} & \colhead{Flag} & \colhead{$N_\mathrm{kin}$\tablenotemark{a}} & \colhead{$N_\mathrm{tot}$\tablenotemark{b}}  & \colhead{$f_\mathrm{w.c.}$\tablenotemark{c}}\\'
printf,lun,'     \colhead{Resolution} & \colhead{} & \colhead{} & \colhead{}  & \colhead{(\%)}'
printf,lun,'  }'
printf,lun,'  \startdata'

kdars = ['N','F','T','O','U','X']
kname = ['Near','Far','Tangent','Outer Galaxy','Unconstrained','Excluded\tablenotemark{d}']

kinfo = constrain.kdist OR constrain.grsmatch OR constrain.knownd OR constrain.parallax
print,fix(total(kinfo))
iki = where(kinfo)
wc = constrain.post
iwc = where(wc)
p = pvec[iki]
c = constrain[iki]
inkd = where(~c.kdist AND ~c.grsmatch,nnkd)
print,nnkd
print,fix(total(wc))
nwc = fix(total(wc))

print,'-------------------------------------------------'
ct=0
kct = fix(total(kinfo))
FOR jj=0,n_elements(kdars)-1 DO BEGIN
   
   fwc = fix(total(kdars[jj] EQ ['N','F','T','O']))
   sfwc = fwc EQ 0 ? '\nodata' : $
          string(total(p.stat.kdar EQ kdars[jj])/nwc*100., format="(F0.1)")
   
   str = '  '+kname[jj]+' & '+kdars[jj]+' & ' + $
         string(fix(total(p.stat.kdar EQ kdars[jj])),format="(I0)") + ' & '+$
         string(fix(total(pvec.stat.kdar EQ kdars[jj])),format="(I0)") + $
         ' & '+$
         sfwc + ' \\'
   printf,lun,str
   ct += fix(total(p.stat.kdar EQ kdars[jj]))
ENDFOR

printf,lun,'  \hline'
printf,lun,'  Total & & '+string(ct,format="(I0)") +' & '+$
       string(n_elements(pvec),format="(I0)") + ' & 100'

printf,lun,'  \enddata'
printf,lun,'  \tablenotetext{a}{Number of objects in the ``kinematic sample''''.}'
printf,lun,'  \tablenotetext{b}{Number of objects from the full Bolocat V2.}'
printf,lun,'  \tablenotetext{c}{Fraction of the well-constrained sources with this KDA resolution.}'
printf,lun,'  \tablenotetext{d}{Object lies in a kinematic avoidance zone (\S\ref{res:kaz}).}'
printf,lun,'\end{deluxetable}'

close,lun
free_lun,lun



;;=================================================================
;;=PLOT=PLOT=PLOT=PLOT=PLOT=PLOT=PLOT=PLOT=PLOT=PLOT=PLOT=PLOT=PLOT
xsize = 3.1
ysize = 3.7
myps,'./grs_paper/glon_distributions.eps',xsize=xsize,ysize=ysize,thick=1.5
!x.omargin = [-0.5,-2.25]
!y.omargin = [-2.5,-3.25]
ygap = 0.035
multiplot,[1,3],/doxaxis,ygap=ygap

cs = 0.55

ineg = where(pvec.glon GE 250)
pvec[ineg].glon -= 360.

wc = constrain.post
iwc = where(wc)
ipc = where(~wc)
p = pvec[ipc]
c = constrain[ipc]

print,'NWC / NPC: ',n_elements(iwc),n_elements(ipc)

cgPlot,[0],[0],/nodata,xr=[222.5,-17.5],yr=[0,950],$
       charsize=cs,xtit='Galactic Longitude  [deg]',$
       ytit='N per 5'+cgSymbol('deg')+' bin'

kazcol = 'blk4'
cgColorfill,[21,!x.crange[1],!x.crange[1],21,21],$
            [0,0,!y.crange[1],!y.crange[1],0],/line_fill,color=kazcol,$
            orientation=45,thick=1
cgColorfill,[70,100,100,70,70],$
            [0,0,500,500,0],/line_fill,color=kazcol,$
            orientation=45,thick=1
cgColorfill,[160,200,200,160,160],$
            [0,0,500,500,0],/line_fill,color=kazcol,$
            orientation=45,thick=1
vline,21,color=kazcol
vline,70,color=kazcol,yr=[0,500]
vline,100,color=kazcol,yr=[0,500]
vline,160,color=kazcol,yr=[0,500]
vline,200,color=kazcol,yr=[0,500]
         

iii = where(p.glon LE 21 OR abs(p.glon - 90) LE 20. OR abs(p.glon - 180) LE 20., niii)
print,'NIII: ',niii
iii = where(p.glon GE 195, niii)
print,'NIII: ',niii

vt = read_mrt('./local/mrt/Shirley_HHT.mrt')

plothist,vt.glon,bin=5,color='red',thick=5,/over


plothist,p.glon,bin=5,color='cyan',thick=2,/over,xarr,yarr
;; cgOplot,xarr,yarr/2.,color='blue',psym=10 
print,!x.crange
p = pvec[iwc]
c = constrain[iwc]
plothist,p.glon,bin=5,/over,color='black',xarr,yarr,thick=2.5
print,max(yarr),max(xarr)


al_legend,/top,/left,box=0,/clear,color=['black','cyan','red'],linestyle=0,$
          linsize=1.5,['Well-Constrained Distances','Unconstrained Distances',$
                      'HHT Obs. (Shirley et al. 2013)'],$
          charsize=1.0*cs,thick=[5,4,10]/2.


cgAxis,xaxis=0,/xst,xtickformat='blank_axis'
cgAxis,xaxis=1,/xst,xtickformat='blank_axis'
cgAxis,yaxis=1,/yst,ytickformat='blank_axis'

multiplot,/doxaxis
;; myps,/done,/mp

;; ;;&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&
;; ;; Another Plot!!!!

;; myps,'./grs_paper/source_histograms.eps',xsize=5
;; !x.omargin = [0,-2]
;; !y.omargin = [-2,-4]
;; multiplot,[1,2],ygap=0.045,/doxaxis

jjj= where(p.glon LT 90, njjj)
print,'NJJJ: ',njjj
p = p[jjj]
c = c[jjj]

;; FAR
ind = where(p.stat.kdar EQ 'F', nfar)
print,nfar
plothist,p[ind].glat,bin=0.1,charsize=cs,xtit='Galacitc Latitude  [deg]',$
         /fill,fcolor='cyan',color='blue',xr=[-1.25,1.25],xarr,yarr,$
         ytit='N per 0.1'+cgSymbol('deg')+' bin',yr=[0,150],ytickinterval=50
print,[min(xarr),max(xarr)]
yf = mpfitpeak(xarr,yarr,A,NTERMS=3)
print,'GLAT, FAR: ',[A,2.355*A[2]]
;; cgOplot,xarr,gauss_1(xarr,A),color='cyan'
bfar = p[ind].glat

;; TAN
ind = where(p.stat.kdar EQ 'T',ntan)
print,ntan
plothist,p[ind].glat,bin=0.1,/over,color='red',xarr,yarr,$
         /fill,/fline,fcolor='red',forientation=45,fthick=1
print,[min(xarr),max(xarr)]
yf = mpfitpeak(xarr,yarr,A,NTERMS=3)
print,'GLAT, TAN: ',[A,2.355*A[2]]
;; cgOplot,xarr,gauss_1(xarr,A),color='red'
btan = p[ind].glat


;; NEAR
ind = where(p.stat.kdar EQ 'N',nnear)
print,nnear
plothist,p[ind].glat,bin=0.1,/over,color='black',xarr,yarr
print,[min(xarr),max(xarr)]
yf = mpfitpeak(xarr,yarr,A,NTERMS=3)
print,'GLAT, NEAR: ',[A,2.355*A[2]]
;; cgOplot,xarr,gauss_1(xarr,A),color='black'
bnear = p[ind].glat

al_legend,/top,/right,box=0,color=['black','red','cyan'],linestyle=0,$
          linsize=1.5,['Near','Tangent','Far'],charsize=1.0*cs

print,!x.crange
kstwo, bfar,btan,D,prob
print,'GLAT, FAR-TAN: ',D,prob
kstwo, bfar,bnear,D,prob
print,'GLAT, FAR-NEAR: ',D,prob
kstwo, bnear,btan,D,prob
print,'GLAT, NEAR-TAN: ',D,prob

cgAxis,xaxis=0,/xst,xtickformat='blank_axis'


multiplot

s = omni_read_cat()
s = s[iwc]
s = s[jjj]

;; FAR
ind = where(p.stat.kdar EQ 'F', nfar)
print,nfar
plothist,alog10(s[ind].flux),bin=0.1,charsize=cs,xtit='Flux Density  [Jy]',$
         /fill,fcolor='cyan',color='blue',xarr,yarr,xst=5,$
         yr=[0,90],ytit='N per 0.1 dex bin'
print,[min(xarr),max(xarr)]
yf = mpfitpeak(xarr,yarr,A,NTERMS=3)
print,'FLUX, FAR: ',[A,10.^(median(alog10(s[ind].flux)))]
;; cgOplot,xarr,gauss_1(xarr,A),color='cyan'
;; vline,median(alog10(s[ind].flux)),linestyle=3,color='cyan'
sfar = alog10(s[ind].flux)

;; TAN
ind = where(p.stat.kdar EQ 'T',ntan)
print,ntan
plothist,alog10(s[ind].flux),bin=0.1,/over,color='red',xarr,yarr,$
         /fill,/fline,fcolor='red',forientation=45,fthick=1
print,[min(xarr),max(xarr)]
yf = mpfitpeak(xarr,yarr,A,NTERMS=3)
print,'FLUX, TAN: ',[A,10.^(median(alog10(s[ind].flux)))]
;; cgOplot,xarr,gauss_1(xarr,A),color='red'
;; vline,median(alog10(s[ind].flux)),linestyle=3,color='red'
stan = alog10(s[ind].flux)

;; NEAR
ind = where(p.stat.kdar EQ 'N',nnear)
print,nnear
plothist,alog10(s[ind].flux),bin=0.1,/over,color='black',xarr,yarr
print,[min(xarr),max(xarr)]
yf = mpfitpeak(xarr,yarr,A,NTERMS=3)
print,'FLUX, NEAR: ',[A,10.^(median(alog10(s[ind].flux)))]
;; cgOplot,xarr,gauss_1(xarr,A),color='black'
;; vline,median(alog10(s[ind].flux)),linestyle=3,color='black'
print,!x.crange
snear = alog10(s[ind].flux)

al_legend,/top,/right,box=0,color=['black','red','cyan'],linestyle=0,$
          linsize=1.5,['Near','Tangent','Far'],charsize=1.0*cs


kstwo, sfar,stan,D,prob
print,'S1.1, FAR-TAN: ',D,prob,(1.d - prob)*100.d
kstwo, sfar,snear,D,prob
print,'S1.1, FAR-NEAR: ',D,prob,(1.d - prob)*100.d
kstwo, snear,stan,D,prob
print,'S1.1, NEAR-TAN: ',D,prob,(1.d - prob)*100.d

;; plothist,alog10(s.flux),bin=0.1,charsize=1.0,xst=5
cgAxis,xaxis=0,/xlog,xr=10.^!x.crange,/save,xtit='Flux Density  [Jy]',$
       charsize=cs,xtickformat='exponent10',/xst
cgAxis,xaxis=1,/xst,/xlog,xtickformat='blank_axis'

myps,/done,/mp


;;(((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((
myps,'./grs_paper/dsun_rgal_distributions.eps',xsize=xsize,ysize=ysize,thick=1.5
;; !x.omargin = [-0,-4]
;; !y.omargin = [-5,-5]
multiplot,[1,3],/doxaxis,ygap=ygap

ikin = where(constrain.kdist OR constrain.grsmatch OR constrain.knownd $
             OR constrain.parallax, nkin)
print,'NKIN: ',nkin
iwc = where(constrain.post)

;; PML Histogram

plothist,pvec[ikin].stat.pml,charsize=cs,xtit='P!dML!n',bin=0.01,yr=[0,425],$
         yminor=4,ytit='N per 0.01 bin'

vline,0.78,color='blk4';,/log
plothist,pvec[iwc].stat.pml,/over,bin=0.01,color='blk4',/fill,$
         fcolor='blk4';,/ylog
plothist,pvec[ikin].stat.pml,/over,bin=0.01

cgAxis,xaxis=0,/xst,xtickformat='blank_axis'
cgAxis,yaxis=1,/yst,ytickformat='blank_axis',yminor=4

al_legend,/top,/right,box=0,linestyle=0,color=['black','blk4'],$
          ['Kinematic Sample','Well-Constrained'],linsize=1.5,charsize=cs*1.0


multiplot,/doxaxis
;;*******************************

BIN = 0.5
stbin = string(BIN,format="(F0.1)")

dsun = pvec[iwc].stat.duse[0]/1.d3

plothist,dsun,bin=bin,charsize=cs,xtit='Heliocentric Distance  [kpc]',$
         ytit='N per '+stbin+' kpc bin'

ind = where(dsun LE 5.5, n5)
print,'DSUN: ',n5,n_elements(dsun),float(n5)/n_elements(dsun)


multiplot,/doxaxis
;;*******************************

;; Compute RGAL from vlsr
restore,'./local/'+conf.survey+'_velocities.sav',/ver

BIN = 0.2
stbin = string(BIN,format="(F0.1)")

vlsr = fltarr(n_elements(v)) - 1000.d
gf = v.grs.flag
igrs = where(gf EQ 1 OR gf EQ 2 OR gf EQ 5 OR gf EQ 6, ngrs)
print,'NGRS:', ngrs
vlsr[igrs] = v[igrs].grs.vlsr
idg  = where(v.vlsr GE -500, ndg)
print,'NDG: ',ndg
vlsr[idg] = v[idg].vlsr

vlsr[where(vlsr LE -500)] = !values.f_nan

print,m4_stat(vlsr[ikin])

;; FOR W-C, compute RGAL from dsun
mw = omni_read_conffile('./conffiles/galactic_params.conf')

rgal = (mw.v0 * sin(pvec.glon*!const.dtor) * mw.r0 / $
        (vlsr + mw.v0*sin(pvec.glon*!const.dtor))) / $
       1.d3 / cos(pvec.glat*!const.dtor)

omni_lbd2rz,pvec[iwc].glon,pvec[iwc].glat,pvec[iwc].stat.duse[0],rgwc,zgwc

help,rgal,rgwc
rgal[iwc] = rgwc / 1.d3

print,m4_stat(rgal,rgal[ikin],rgwc/1.d3)


plothist,rgal[ikin],bin=bin,charsize=cs,xtit='Galactocentric Radius  [kpc]',$
         ytit='N per '+stbin+' kpc bin',yminor=4

;; Make filled plot for WC sample only
plothist,rgal[iwc],bin=bin,/over,color='blk4',/fill,fcolor='blk4'
plothist,rgal[ikin],bin=bin,/over

cgAxis,xaxis=0,/xst,xtickformat='blank_axis'

al_legend,/top,/right,box=0,linestyle=0,color=['black','blk4'],$
          ['Kinematic Sample','Well-Constrained'],linsize=1.5,charsize=cs*1.0



myps,/done,/mp

;;&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&
;; Yet Another Plot!!!!
myps,'./grs_paper/flux_distance.eps',xsize=5

p = pvec[iwc]
c = constrain[iwc]
s = omni_read_cat()
s = s[iwc]



cgPlot,s.flux,p.stat.duse[0]/1.d3,psym=16,symsize=0.4,/xlog,/xst,$
       xtickformat='exponent10',xtit='Flux Density  [Jy]',$
       ytit='Heliocentric Distance  [kpc]',charsize=cs,xr=10.^[ -1.15,2.25]

hist = hist_2d(alog10(s.flux),p.stat.duse[0]/1.d3,$
               bin1=0.1,min1=-1.15,max1=2.25,$
               bin2=0.5,min2=0,max2=20)
;; help,hist
;; cgLoadct,13
;; pr = set_plot_range(hist)
;; print,'PLOT RANGE: ',pr
;; plotimage,hist,range=pr,xst=4,yst=4
;; cgAxis,yaxis=0,/save,yr=[0,20],ytit='Heliocentric Distance  [kpc]',charsize=1.0
;; cgAxis,yaxis=1,ytickformat='blank_axis',/yst

;; cgColorbar,/vert,/right,range=pr,ytickinterval=2

print,'2D:'
print,!x.crange,!y.crange


myps,/done
;;=PLOT=PLOT=PLOT=PLOT=PLOT=PLOT=PLOT=PLOT=PLOT=PLOT=PLOT=PLOT=PLOT
;;;;=================================================================



















;;====================================
;; Table Header
openw,lun,'./grs_paper/dpdf_matrix.tex',/get_lun

printf,lun,'%%=================================%%'
printf,lun,'%%                                 %%'
printf,lun,'%%   Table: DPDF Matrix            %%'
printf,lun,'%%   Label: table:matrix           %%'
printf,lun,'%%                                 %%'
printf,lun,'%%=================================%%'
printf,lun,'\begin{deluxetable}{'+strjoin(replicate('c',ndpdf+2))+'}'
printf,lun,'  \tablecolumns{'+string(ndpdf+2,format="(I0)")+'}'
printf,lun,'  \tablewidth{0pt}'
printf,lun,'  \tabletypesize{\footnotesize}'
printf,lun,'  \tablecaption{DPDF Overlap Matrix\label{table:matrix}}'
printf,lun,'  \tablehead{'
printf,lun,'    \colhead{} & \multicolumn{2}{c}{ $\mathcal{L}(v_\mathrm{LSR},l,b;d_\sun)$\tablenotemark{a}} & \colhead{} & \multicolumn{5}{c}{Prior DPDFs} & \colhead{} \\'
printf,lun,'    \cline{2-3} \cline{5-9}'
str = '    \colhead{} & '
FOR ii=0,ndpdf-1 DO BEGIN
   IF ii GT 0 THEN str += ' & '
   str += '\colhead{'+translate_dpdf_tag(dnames[ii])+'}'
   IF ii EQ 1 THEN BEGIN
      print,translate_dpdf_tag(dnames[ii])
      str += ' & '
   ENDIF
ENDFOR
printf,lun,str
printf,lun,'  }'
printf,lun,'  \startdata'


matrix = dblarr(ndpdf,ndpdf)
postma = dblarr(ndpdf,ndpdf)

;; Make Matrix:
FOR ii=0,ndpdf-1 DO BEGIN
   str = '  '+translate_dpdf_tag(dnames[ii])
   IF ii GT 0 THEN str += strjoin(replicate(' & ',ii))
   IF ii GT 1 THEN str += ' & '
   FOR jj=ii, ndpdf-1 DO BEGIN
      str += ' & '
      matrix[ii,jj] = total(constrain.(ii+1) AND constrain.(jj+1))
      postma[ii,jj] = total(constrain.(ii+1) AND constrain.(jj+1) AND $
                            constrain.post)
      str += string(matrix[ii,jj],format="(I0)")
      IF jj EQ 1 THEN str += ' & '
   ENDFOR
   str += ' \\'
   printf,lun,str
ENDFOR



printf,lun,'  \enddata'
printf,lun,'  \tablenotetext{a}{Kinematic distance likelihood; dense gas and \thco\ are mutually exclusive.}'
printf,lun,'\end{deluxetable}'

close,lun
free_lun,lun



;;======================================================================
;; Prior DPDF overlap / well-constrained fraction


;;====================================
;; Table Header
openw,lun,'./grs_paper/priors_overlap.tex',/get_lun

printf,lun,'%%=================================%%'
printf,lun,'%%                                 %%'
printf,lun,'%%   Table: Prior DPDF Overlap     %%'
printf,lun,'%%   Label: table:priors           %%'
printf,lun,'%%                                 %%'
printf,lun,'%%=================================%%'
printf,lun,'\begin{deluxetable*}{l'+strjoin(replicate('c',ndpdf-3))+'}'
printf,lun,'  \tablecolumns{'+string(ndpdf-2,format="(I0)")+'}'
printf,lun,'  \tablewidth{0pt}'
printf,lun,'  \tabletypesize{\footnotesize}'
printf,lun,'  \tablecaption{Well-Constrained Fraction for Overlapping Priors\label{table:priors}}'
printf,lun,'  \tablehead{'
str = '    \colhead{} & '
FOR ii=2,ndpdf-2 DO BEGIN
   str += '\colhead{'+translate_dpdf_tag(dnames[ii])+'}'
   IF ii NE ndpdf-2 THEN str += ' & '
ENDFOR
printf,lun,str
printf,lun,'  }'
printf,lun,'  \startdata'



;; Make Matrix:
FOR ii=2,ndpdf-2 DO BEGIN
   str = '  '+translate_dpdf_tag(dnames[ii])
   IF ii GT 2 THEN str += strjoin(replicate(' & ',ii-2))
   FOR jj=ii, ndpdf-2 DO BEGIN
      str += ' & '
      val = postma[ii,jj] / matrix[ii,jj]
      str += finite(val) ? string(val,format="(F0.3)") : '\nodata'
      IF jj EQ 1 THEN str += ' & '
   ENDFOR
   IF ii NE ndpdf-2 THEN str += ' \\'
   printf,lun,str
ENDFOR



printf,lun,'  \enddata'
printf,lun,'\end{deluxetable*}'

close,lun
free_lun,lun



;;==========================================
;; Another Table!

openw,lun,'./grs_paper/dpdf_applicability.tex',/get_lun

printf,lun,'%%=================================%%'
printf,lun,'%%                                 %%'
printf,lun,'%%   Table: DPDF Applicability     %%'
printf,lun,'%%   Label: table:dpdfs            %%'
printf,lun,'%%                                 %%'
printf,lun,'%%=================================%%'
printf,lun,'\begin{deluxetable*}{lrrccc}'
printf,lun,'  \tablecolumns{6}'
printf,lun,'  \tablewidth{0pt}'
printf,lun,'  \tabletypesize{\footnotesize}'
printf,lun,'  \tablecaption{Application of DPDFs to BGPS V2 Sources\label{table:dpdfs}}'
printf,lun,'  \tablehead{'
printf,lun,'    \colhead{} & \colhead{} & \colhead{} & \multicolumn{3}{c}{Well-Constrained Sources Only} \\'
printf,lun,'    \cline{4-6}'
printf,lun,'    \colhead{DPDF}  & \colhead{$N$} & \colhead{$N_\mathrm{wc}$} & \colhead{Fraction of} & \colhead{Fraction of} & \colhead{Fraction of} \\'
printf,lun,'    \colhead{} & \colhead{} & \colhead{} & \colhead{Method (\%)} & \colhead{w.c.\tablenotemark{a} (\%)} & \colhead{Full V2.1 (\%)}'
printf,lun,'  }'
printf,lun,'  \startdata'


wckd = constrain.(1) AND constrain.post

FOR ii=0,ndpdf-1 DO BEGIN
   
   IF ii EQ ndpdf-1 THEN printf,lun,'\hline'
   
   wc = constrain.(ii+1) AND constrain.post
   
   str = '  '+translate_dpdf_tag(dnames[ii])+ ' & '+$
         string(total(constrain.(ii+1)),total(wc),format="(I0,' & ',I0)")
   str += ' & '
   CASE 1 OF
      ii EQ ndpdf-1: this = '\nodata & '
      total(wc) / total(constrain.(ii+1)) EQ 1: this = '100 & '
      ELSE: this = string(total(wc) / $
                          total(constrain.(ii+1))*100.,format="(F0.1)")+' & '
   ENDCASE
   str += this
   CASE 1 OF
      ;;ii EQ ndpdf-1: this = '\nodata'
      total(wc) / total(constrain.post) EQ 1: this = '100'
      ELSE: this = string(total(wc) / $
                          total(constrain.post)*100.,format="(F0.1)")
   ENDCASE
   str += this
   str += ' & '+string(total(wc) / n_elements(pvec)*100.,format="(F0.1)")
   IF ii NE ndpdf-1 THEN str += ' \\'
   printf,lun,str
   
ENDFOR

printf,lun,'  \enddata'
printf,lun,'  \tablenotetext{a}{Fraction of the set of well-constrained sources.}'
printf,lun,'\end{deluxetable*}'


close,lun
free_lun,lun





END
